\documentclass[10pt,letterpaper]{article}

\usepackage[latin1]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{moreverb}
\usepackage{graphicx}
\usepackage{color}
\usepackage{fancyheadings}
\usepackage{listings}
\usepackage{algorithmic}
\usepackage{algorithm}

\graphicspath{{../figures/}}

%-----------------------------------------------------------------------
% Text dimensions 
%-----------------------------------------------------------------------
\setlength{\textheight}{22cm}
\addtolength{\textwidth}{56pt} % Aproxx. ~16cm
%
% Becuase I don't need marginal notes, then: 
%
\setlength{\marginparsep}{0pt}
\setlength{\marginparwidth}{0pt}
%
% Header is equal to \textwdith
% This is used when \marginparsep and \marginparwith are not zero
%
%\addtolength{\headwidth}{\marginparsep}
%\addtolength{\headwidth}{\marginparwidth}
%
% Header is equal to \textwdith
%
\setlength{\headwidth}{\textwidth}
%
% We have to substrac some points on even pages,
% because we add some points to \textwidth
%
\addtolength{\evensidemargin}{-56pt}
%\addtolength{\oddsidemargin}{-10pt}

\setlength{\voffset}{-1cm}


%-----------------------------------------------------------------------
% Style of page and headings
%-----------------------------------------------------------------------
\pagestyle{fancyplain}
\renewcommand{\sectionmark}[1]
	{\markright{\thesection\ #1}}
\lhead[\fancyplain{}{\bf\thepage}]
	{\fancyplain{}{\bfseries\rightmark}}
\rhead[\fancyplain{}{\bfseries\leftmark}]
	{\fancyplain{}{\bf\thepage}}
\cfoot{}
%\lfoot[\fancyplain{}{\bfseries\thepage}]{\fancyplain{}{}}
%\rfoot[\fancyplain{}{}]{\fancyplain{}{\bfseries\thepage}}

\newcommand{\mathbi}[1]{\textbf{\em #1}}
\newcommand{\mbf}[1]{\mathbf{#1}}
\newcommand{\mbi}[1]{\textbf{\em #1}}
\newcommand{\mbiu}[1]{{_\textbf{\em #1}}}
\newcommand{\mcons}[1]{\mbox{\textsf{#1}}}
\newcommand{\tensor}[1]{\underline{\underline{#1}}}


\lstset{language=c++}
\lstset{commentstyle=\textit}


\author{Luis M. de la Cruz}

\title{TUNA: Meshes and Storage}

\begin{document}

\maketitle

\tableofcontents

\section{Introduction}

\subsection{}

\section{Meshes}

The class \textsf{Mesh} is an abstract class that is the root of the hierarchy for all 
kind of meshes. This class store the nodes, but not the conectivity, which depends
on the type of the mesh to be constructed, so only specialized classes can construct it. 
Curiosly Recurring Template Pattern is used to delegate responsabilities to more specialized 
meshes, for example to calculate coordinates of the nodes and conectivity. By now, only 
structured meshes in Cartesian coordinates are implemented, but more complex and unstructured 
meshes should be possible to add, with a little bit of effort.

\begin{figure}[h!]
$$\includegraphics[width=14cm]{meshes.png}$$
\caption{Mesh representation in TUNA. The FVM requires the nodes (a)
and the volumes (b). In TUNA the number of volumes (\textsf{Nvx} and \textsf{Nvy}) 
includes the boundaries.}
\label{fig:meshes}
\end{figure}

\subsection{Structured meshes}

The structured meshes in TUNA are represented as shown in figure \ref{fig:meshes}. In the Finite
Volume Method (FVM) nodes and volumes are required. In the example of the figure \ref{fig:meshes}, 
the number of nodes are \textsf{Nx} and \textsf{Ny} in $x$ and $y$ axis respectively. The number of
actual control volumes is $(\mcons{Nx} - 1) \times (\mcons{Ny} - 1)$, see figure \ref{fig:meshes}(b).
Because we need also to account for the boundaries (white points in figure \ref{fig:meshes}(b)), 
the mesh objects will always calculate the number of volumes as follows:
$\mcons{Nvx} = (\mcons{Nx} - 1) + 2 = \mcons{Nx} + 1$, which represents the number of 
{\em "volumes"} in the x-axis. The same is done in the y-axis. The class \textsf{StructuredMesh}
has many member functions that calculate and returns important data of the mesh objects.

\subsubsection{Uniform meshes}

The simplest case of structured meshes is when everything is uniform. There is an adaptor called
\textsf{Uniform}. This is used only to declare and define uniform structured meshes. Next example
construct the object \textsf{mymesh} which represents a mesh of this type:

\strut

\begin{scriptsize}
\begin{lstlisting}[frame=trbl]{}
int Nx, Ny; 
double Lx, Ly; 
StructuredMesh <Uniform<double, 2 > > mymesh (Lx, Nx, Ly, Ny); 
\end{lstlisting}
\end{scriptsize}

\strut

Here \textsf{Nx} and \textsf{Ny} are the number of nodes,
and \textsf{Lx} and \textsf{Ly} are the lenghts of the domain in $x$ and $y$-axis respectively.
These data are user defined.
The \textsf{Uniform} adaptor can be parametrized with the precision (\textsf{float}, 
\textsf{double} or \textsf{long double}) and with the dimension (1, 2 or 3).

\subsection{Non--uniform meshes}

There is also the possibility to construct non--uniform meshes using the \textsf{NonUniform}
adaptor. This adaptor implements several functions to calculate the distribution of points
along the three axis. What is required from the user is the definition of a function to 
distribute the points. This is done using {\em functor}'s:

\strut

\begin{scriptsize}
\begin{lstlisting}[frame=trbl]{}
// Functor to define the point distribution.
class MyFunction {      
  double lambda; // distribution factor 
public:                                    
  MyFunction(double l) : lambda(l) { }   
  
  double operator() (double epsi) {      
    return ( (exp(lambda * epsi) - 1.0) / 
	     (exp(lambda       ) - 1.0) ); 
  }
};
\end{lstlisting}
\end{scriptsize}

Then, to create a non--uniform mesh object, we need to do the following:

\begin{scriptsize}
\begin{lstlisting}[frame=trbl]{}
StructuredMesh<NonUniform<double, 2> > mesh(Lx, Nx,Ly, Ny);
\end{lstlisting}
\end{scriptsize}                                  

Finally, the point distribution is calculated calling the \textsf{calcPoints()} 
member function, passing the \textsf{MyFunction} functor as argument:  
\begin{scriptsize}
\begin{lstlisting}[frame=trbl]{}
mesh.calcPoints(X, MyFunction(lambda) , 0       , Lx * 0.5, Nx / 2 + 1);            
mesh.calcPoints(X, MyFunction(-lambda), Lx * 0.5, Lx      , Nx / 2 + 1); 
mesh.calcPoints(Y, MyFunction(lambda) , 0       , Ly * 0.5, Ny / 2 + 1);
mesh.calcPoints(Y, MyFunction(-lambda), Ly * 0.5, Ly      , Ny / 2 + 1);       
\end{lstlisting}
\end{scriptsize}

\strut

The arguments are: the axis in which the point distribution will be calculated,
the functor, the starting and finish coordinate where the distribution will
be calculated, and the number of points to be distributed. The result of the 
previous example is shown in figure \ref{fig:meshesNonuniform}, where you can
see different distribution in each half of the sides of the domain.

\begin{figure}[h!]
$$\includegraphics[width=7cm]{meshesNonuniform.png}$$
\caption{Non--uniform mesh with several distribution of points.}
\label{fig:meshesNonuniform}
\end{figure}

\section{Storage}



The equations will be discretized using a mesh of volumes, see figure \ref{fig:meshes}(b). For example, the
discretized Laplace equation yields a matrix like:
%
%\begin{figure}[h!]
%$$\includegraphics[width=6cm]{meshFVM.png}$$
%\caption{Mesh for the discretization using FVM. The indexes $(i,j)$ represent a
%volume on the mesh.}
%\label{fig:meshFVM}
%\end{figure}

\begin{scriptsize}
\begin{equation}
\begin{array}{c|c|cccccc|cccccc|cccccc|cccccc|}
i  &   & 1 & 2 & 3 & 4 & 5 & 6 & 1 & 2 & 3 & 4 & 5 & 6 & 1 & 2 & 3 & 4 & 5 & 6 & 1 & 2 & 3 & 4 & 5 & 6 \\
\hline
   & j & \multicolumn{6}{c|}{1} & \multicolumn{6}{c|}{2} & \multicolumn{6}{c|}{3} & \multicolumn{6}{c|}{4} \\
\hline
1 & 1 & a & b &   &   &   &   & d &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   &  \\
2 & 1 & c & a & b &   &   &   &   & d &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   &  \\
3 & 1 &   & c & a & b &   &   &   &   & d &   &   &   &   &   &   &   &   &   &   &   &   &   &   &  \\
4 & 1 &   &   & c & a & b &   &   &   &   & d &   &   &   &   &   &   &   &   &   &   &   &   &   &  \\
5 & 1 &   &   &   & c & a & b &   &   &   &   & d &   &   &   &   &   &   &   &   &   &   &   &   &  \\
6 & 1 &   &   &   &   & c & a &   &   &   &   &   & d &   &   &   &   &   &   &   &   &   &   &   &  \\
\hline
1 & 2 & e &   &   &   &   &   & a & b &   &   &   &   & d &   &   &   &   &   &   &   &   &   &   &  \\
2 & 2 &   & e &   &   &   &   & c & a & b &   &   &   &   & d &   &   &   &   &   &   &   &   &   &  \\
3 & 2 &   &   & e &   &   &   &   & c & a & b &   &   &   &   & d &   &   &   &   &   &   &   &   &  \\
4 & 2 &   &   &   & e &   &   &   &   & c & a & b &   &   &   &   & d &   &   &   &   &   &   &   &  \\
5 & 2 &   &   &   &   & e &   &   &   &   & c & a & b &   &   &   &   & d &   &   &   &   &   &   &  \\
6 & 2 &   &   &   &   &   & e &   &   &   &   & c & a &   &   &   &   &   & d &   &   &   &   &   &  \\
\hline
1 & 3 &   &   &   &   &   &   & e &   &   &   &   &   & a & b &   &   &   &   & d &   &   &   &   &   \\
2 & 3 &   &   &   &   &   &   &   & e &   &   &   &   & c & a & b &   &   &   &   & d &   &   &   &   \\
3 & 3 &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a & b &   &   &   &   & d &   &   &   \\
4 & 3 &   &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a & b &   &   &   &   & d &   &   \\
5 & 3 &   &   &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a & b &   &   &   &   & d &   \\
6 & 3 &   &   &   &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a &   &   &   &   &   & d \\
\hline
1 & 4 &   &   &   &   &   &   &   &   &   &   &   &   & e &   &   &   &   &   & a & b &   &   &   &   \\
2 & 4 &   &   &   &   &   &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a & b &   &   &   \\
3 & 4 &   &   &   &   &   &   &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a & b &   &   \\
4 & 4 &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a & b &   \\
5 & 4 &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a & b \\
6 & 4 &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   &   & e &   &   &   &   & c & a \\
\hline
\end{array}
\label{eq:laplace}
\end{equation}
\end{scriptsize}

In the above matrix, the $i$ and $j$ indexes identify the volumes on the mesh, figure \ref{fig:meshes}(b).
The matrices resulting from the discretization are banded, and for the 2D case we get 5 diagonals. 
For 1D and 3D the matrices have 3 and 7 diagonals respectively. This kind of matrices can be stored 
using an array for each diagonal. In the example, the arrays are as follows:

\begin{scriptsize}
\begin{equation}
\begin{array}{rr}
i  &
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
 1 & 2 & 3 & 4 & 5 & 6 &  1 & 2 & 3 & 4 & 5 & 6 & 1 & 2 & 3 & 4 & 5 & 6 & 1 & 2 & 3 & 4 & 5 & 6 \\
\end{array} \\
\hline
j  &
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
 1 & 1 & 1 & 1 & 1 & 1 &  2 & 2 & 2 & 2 & 2 & 2 & 3 & 3 & 3 & 3 & 3 & 3 & 4 & 4 & 4 & 4 & 4 & 4 \\
\hline
\end{array} \\
\hline
a  \rightarrow &
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
\hline
\end{array} \\
b \rightarrow &
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
 1 & 1 & 1 & 1 & 1 & \textcolor{red}{0} &
 1 & 1 & 1 & 1 & 1 & \textcolor{red}{0} &
 1 & 1 & 1 & 1 & 1 & \textcolor{red}{0} &
 1 & 1 & 1 & 1 & 1 & \textcolor{red}{0}\\
\hline
\end{array} \\
c \rightarrow &
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
 \textcolor{red}{0} & 1 & 1 & 1 & 1 & 1 & 
 \textcolor{red}{0} & 1 & 1 & 1 & 1 & 1 & 
 \textcolor{red}{0} & 1 & 1 & 1 & 1 & 1 & 
 \textcolor{red}{0} & 1 & 1 & 1 & 1 & 1  \\
\hline
\end{array} \\
d \rightarrow &
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 &
 \textcolor{red}{0} & \textcolor{red}{0} & \textcolor{red}{0} & \textcolor{red}{0} & \textcolor{red}{0} & \textcolor{red}{0} \\
\hline
\end{array} \\
e \rightarrow &
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
 \textcolor{red}{0} & \textcolor{red}{0} & \textcolor{red}{0} & \textcolor{red}{0} & \textcolor{red}{0} & \textcolor{red}{0} &  
 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
\hline
\end{array} \\
\end{array}
\label{eq:storagePenta}
\end{equation}
\end{scriptsize}

In this format, the number 1 in the arrays $a,b,c,d$ and $e$ means that a value other 
than zero is stored and correspond to the value of the coefficient calculated in the 
node $(i,j)$ of the mesh. Using this strategy, the storage is reduced from
$(\mcons{Nx} - 1)^2 \times (\mcons{Ny} - 1)^2$ to 
$5 * (\mcons{Nx} - 1) \times (\mcons{Ny} - 1)$. 
Defining \textsf{N} as the total number of real volumes in the mesh, i.e. 
$\mcons{N} = (\mcons{Nx} - 1) * (\mcons{Ny} - 1)$, then the order is reduced 
from $\mathcal{O}(\mcons{N}^2)$ to $\mathcal{O}(\mcons{N})$.
This storage is useful and efficient for the calculation of FVM coefficients in combination with
the TDMA \cite{} solver. For other kind of solvers it is required to change the storage format
in order to take advantage of the sparsity of the matrix and the operations needed for each 
solver (e.g. matrix-vector product).
Using this format we waste some memory: for arrays $b$ and $c$
we miss $\mcons{Ny} - 1$ places; and for arrays $d$ and $e$ we miss $\mcons{Nx} - 1$ places. In total, we
are wasting $2 * ( \mcons{Nx} + \mcons{Ny} - 2)$ places in the storage.

The format (\ref{eq:storagePenta}) is implemented in the class \textsf{DiagonalMatrix} and is
intended to be used in conjunction with the \textsf{TDMA()} solver. The next table resume this 
storage format.

\begin{table}[h!]
\begin{center}
\begin{scriptsize}
\begin{tabular}{|c|c|c|c|}
\hline
 & Class declaration & Total space & Wasted space \\
\hline
1D & \verb|DiagonalMatrix < prec, 1>| & 3 * (\textsf{Nx} - 1) &  2 \\
\hline
2D & \verb|DiagonalMatrix < prec, 2>|  & 5 * (\textsf{Nx} - 1) $\times$ (\textsf{Ny} - 1) &
2 * ( \textsf{Nx + Ny} - 2) \\
\hline
3D & \verb|DiagonalMatrix < prec, 3>| & 7 * (\textsf{Nx} - 1) $\times$ (\textsf{Ny} - 1) $\times$ (\textsf{Nz} - 1) & 
2 * (\textsf{Nx} - 1)(\textsf{Nz} - 1) + \\
& & & 2 * (\textsf{Ny} - 1)(\textsf{Nz} - 1) + \\
& & & 2 *(\textsf{Nx} - 1)(\textsf{Ny} - 1) )\\
\hline
\end{tabular}
\end{scriptsize}
\end{center}
\caption{\textsf{DiagonalMatrix storage}. \textsf{Nx}, \textsf{Ny} and \textsf{Nz} are the number of nodes
in $x$, $y$, and $z$ axis respectively. \textsf{prec} can be \textsf{float}, \textsf{double} or 
\texttt{long double}.} \label{tab:diagonalMatrix}
\end{table}

The wasted space percentage decrease as the number of volumes grows. For example, if 
$\mcons{Nx}-1 = \mcons{Ny}-1 = \mcons{Nz}-1 = \mcons{n}$, then this percentage is of the
order $\mathcal{O}(1/\mcons{n})$.
The implementation of \textsf{DiagonalMatrix} class for the two dimensional case is as follows:

\strut

\begin{scriptsize}
\begin{lstlisting}[frame=trbl]{}
template<typename Tprec>
class DiagonalMatrix<Tprec, 2> {
public:  
/*! Arrays mapped to a mesh in 2D */
  typedef typename TunaArray<Tprec, 2>::huge ScalarField;
  ScalarField a, b, c, d, e;
  
  DiagonalMatrix(int ni, int nj) {
    a.resize(ni, nj);
    b.resize(ni, nj);
    c.resize(ni, nj);
    d.resize(ni, nj);
    e.resize(ni, nj);
  }
  ~DiagonalMatrix() { }
};
\end{lstlisting}
\end{scriptsize}

\strut

This implementation is very simple, it only declares five two dimensional \textsf{TunaArray}'s,
which is only a wrapper of the blitz++ arrays \cite{}. Observe that in this implementation
the arrays (\textsf{a, b, c, d}, and \textsf{e}), are public. This is done just for efficiency. 
The parameters required for this implementation are the precision and the dimension.

A declaration of a \textsf{DiagonalMatrix} object would be as follows:

\strut

\begin{scriptsize}
\begin{lstlisting}[frame=trbl]{}
DiagonalMatrix <double, 2> A(Nx - 1, Ny - 1); // Nx number of nodes in x-axis 
                                              // Ny number of nodes in y-axis
\end{lstlisting}
\end{scriptsize}

\strut

This declaration constructs five two dimensional arrays of size $(\mcons{Nx} - 1) \times (\mcons{Ny} - 1)$.
%In order to fill the matrix (\ref{eq:laplace}) using the format ( \ref{eq:storagePenta}) we just can do
%something like:
%
%\strut
%
%\begin{scriptsize}
%\begin{lstlisting}[frame=trbl]{}
%for(int i = 0; i < Nx; i++)
%   for(int j = 0; i < Ny; j++) {
%     A.a(i,j) = ... // An FVM approximation.
%                ... // The same for the other arrays.
%   }
%\end{lstlisting}
%\end{scriptsize}
%
%\strut
%
%This is very convenient, because now we have a direct correspondence between the matrix and the
%mesh: e.g. to modify the coefficient \textsf{b} corresponding to the node \textsf{(i,j)} of the mesh 
%we just do \textsf{b(i,j) = somevalue}.

\section{General Equation}

The FVM discretizaci\'on of a PDE in three dimension gives as a result a system like 

\begin{equation}
\mbi{a}_\mbiu{P} \psi_\mbiu{P} = \mbi{a}_\mbiu{E} \psi_\mbiu{E} + \mbi{a}_\mbiu{W} \psi_\mbiu{W} +
\mbi{a}_\mbiu{N} \psi_\mbiu{N} + \mbi{a}_\mbiu{S} \psi_\mbiu{S} +
\mbi{a}_\mbiu{F} \psi_\mbiu{F} + \mbi{a}_\mbiu{B} \psi_\mbiu{B} + \mbi{s}_\mbiu{P}
\label{eq:FVM01}
\end{equation}

\noindent
where the shape of coefficients $\mbi{a}$'s and $\mbi{s}_\mbiu{P}$ depends on the numerical scheme used 
to approximate each term of the PDE \cite{}. 
In the \textsf{GeneralEquation} class, these FVM coefficients are defined as arrays with names
\textsf{aP}, \textsf{aE}, \textsf{aW}, \textsf{aN}, \textsf{aS}, \textsf{aF}, \textsf{aB} and \textsf{sp}. 
These arrays are not constructed, but make references to the arrays defined in \textsf{DiagonalMatrix}.
In two dimensions, see (\ref{eq:laplace}), we have

\begin{center}
\begin{tabular}{ccccc}
$\mbi{a}_\mbiu{P} \rightarrow$ \textsf{a} & $\mbi{a}_\mbiu{E} \rightarrow$ \textsf{b} & 
$\mbi{a}_\mbiu{W} \rightarrow$ \textsf{c} & $\mbi{a}_\mbiu{N} \rightarrow$ \textsf{d} & 
$\mbi{a}_\mbiu{S} \rightarrow$ \textsf{e} 
\end{tabular}
\end{center}

These associations means that, for instance, when we modify $\mbi{a}_\mbiu{P}$, we are 
actually modifying the array \textsf{DiagonalMatrix::a}, and so on. 

These associations give us an easy way to relate the volumes numeration of the mesh with
the arrays indexes. For example, figure \ref{fig:malla1d} shows a simple mesh of six internal
volumes, where the equation for the volume 3 is displayed with the contribution from its neighbor
volumes.

\begin{figure}[h!]
$$\includegraphics[width=7cm]{malla1d.png}$$
\caption{Relation between the mesh and the arrays.}
\label{fig:malla1d}
\end{figure}

In figure \ref{fig:malla1d} the arrays $\mbi{a}$'s and $\mbi{s}_\mbiu{P}$ require the value stored
in the entry 3, which is the index of the volume we are working on. On the other hand,
the entries required for the array $\psi$, the unknown, are those at the \textit{west}
and \textit{east} of $\mbi{P}$, i.e. entries 2 and 4, respectively.
Then, the linear system for the mesh of the figure \ref{fig:malla1d} is written as follows

{\footnotesize{
\begin{equation}
\left(
\begin{array}{cccccc}
\mbi{a}_\mbiu{P}(1) & -\mbi{a}_\mbiu{E}(1) & 0 & 0 & 0 & 0 \\
-\mbi{a}_\mbiu{W}(2) & \mbi{a}_\mbiu{P}(2) & -\mbi{a}_\mbiu{E}(2) & 0 & 0 & 0 \\
0 & -\mbi{a}_\mbiu{W}(3) & \mbi{a}_\mbiu{P}(3) & -\mbi{a}_\mbiu{E}(3) & 0 & 0 \\
0 & 0 & -\mbi{a}_\mbiu{W}(4) & \mbi{a}_\mbiu{P}(4) & -\mbi{a}_\mbiu{E}(4) & 0  \\
0 & 0 & 0 & -\mbi{a}_\mbiu{W}(5) & \mbi{a}_\mbiu{P}(5) & -\mbi{a}_\mbiu{E}(5)  \\
0 & 0 & 0 & 0 & -\mbi{a}_\mbiu{W}(6) & \mbi{a}_\mbiu{P}(6) \\
\end{array}\right) \left(
\begin{array}{c}
\psi(1) \\
\psi(2) \\
\psi(3) \\
\psi(4) \\
\psi(5) \\
\psi(6) \\
\end{array} \right) = \left(
\begin{array}{c}
\mbi{s}_\mbiu{P}(1) \\
\mbi{s}_\mbiu{P}(2) \\
\mbi{s}_\mbiu{P}(3) \\
\mbi{s}_\mbiu{P}(4) \\
\mbi{s}_\mbiu{P}(5) \\
\mbi{s}_\mbiu{P}(6) \\
\end{array} \right)
\label{eq:mat2d}
\end{equation}
}}

As in the two dimensional case, here we make the following associations:

\begin{center}
\begin{tabular}{ccc}
$\mbi{a}_\mbiu{P} \rightarrow$ \textsf{a} & $\mbi{a}_\mbiu{E} \rightarrow$ \textsf{b} & 
$\mbi{a}_\mbiu{W} \rightarrow$ \textsf{c} 
\end{tabular}
\end{center}

Note that there is a minus sign in the $\mbi{a}_\mbiu{W}$ and $\mbi{a}_\mbiu{E}$ coefficients
in the system (\ref{eq:mat2d}). This fact must be taken into account when the system be solved.
For example, the tridiagonal system (\ref{eq:mat2d}) can be solved using the Thomas
algorithm, which using the association above mentioned and taken in to account the sign of
$\mbi{a}_\mbiu{W}$ and $\mbi{a}_\mbiu{E}$, is written as shown in algorithm \ref{alg:thomas}.

\begin{algorithm}[h!]
\caption{Thomas}\label{alg:thomas}
{\footnotesize{
\begin{algorithmic}[1]
  \STATE INPUT: $\mbi{a}_\mbiu{P}$, $\mbi{a}_\mbiu{E}$, $\mbi{a}_\mbiu{W}$, $\mbi{s}_\mbiu{P}$
  \STATE OUTPUT: Solution $\psi$
  \STATE $a \leftarrow \mbi{a}_\mbiu{P}$ ; 
  $b \leftarrow \mbi{a}_\mbiu{E}$;
  $c \leftarrow \mbi{a}_\mbiu{W}$; 
  $d \leftarrow \mbi{s}_\mbiu{P}$
  \STATE $\displaystyle P_1 \leftarrow \frac{b(1)}{a(1)}$ 
  \STATE $\displaystyle Q_1 \leftarrow \frac{d(1)}{a(1)}$
  
  \FOR{$k$ = 2 to $N$}     
    \STATE $\displaystyle P_k \leftarrow 
    \frac{b(k)}{a(k) - c(k) P_{k-1}}$
    \STATE $\displaystyle Q_k \leftarrow 
    \frac{d(k) + c(k) Q_{k-1}}{a(k) - c(k) P_{k-1}} $
     (The $+$ in the numerator is because the $\mbi{a}_\mbiu{E}$ and $\mbi{a}_\mbiu{W}$ have a $-$ in the system)
  \ENDFOR
  \STATE $\psi_N \leftarrow Q_N $
 
  \FOR{$k = N-1$ to 1} %(Backward substitution)
    \STATE $\psi_k \leftarrow Q_k + P_k \psi_{k+1}$ (The $+$ is because $P_k$ carry a $-$ from $\mbi{a}_\mbiu{E}$)
  \ENDFOR
\end{algorithmic}
}}
\end{algorithm}


In two dimensions, the FVM discrete equations are written as follows

{\small{
\begin{eqnarray}
\mbi{a}_\mbiu{P}(i,j) \psi(i,j) & = & \mbi{a}_\mbiu{E}(i,j) \psi(i+1,j) + 
                                      \mbi{a}_\mbiu{W}(i,j) \psi(i-1,j) + \label{eq:FVM02}\\ 
                                &   & \mbi{a}_\mbiu{N}(i,j) \psi(i,j+1) + 
                                      \mbi{a}_\mbiu{S}(i,j) \psi(i,j-1) +  
                                      \mbi{s}_\mbiu{P}(i,j) \nonumber 
\end{eqnarray}
}}

To solve this system it is possible to sweep the mesh shown in figure \ref{fig:meshes} by lines in $x$ axis or
in $y$ axis. Because the arrays are stored by rows in memory (the C++ way), it is most efficient to do this swept 
in the $y$ direction. The system (\ref{eq:FVM02}) is transformed into a tridiagonal one considering $\psi(i+1,j)$
and $\psi(i-1,j)$ as known, in such a way we can write

{\small{
\begin{eqnarray}
\mbi{a}_\mbiu{P}(i,j) \psi^{l+1}(i,j) - \mbi{a}_\mbiu{N}(i,j) \psi^{l+1}(i,j+1) - \mbi{a}_\mbiu{S}(i,j) \psi^{l+1}(i,j-1)
& = & \mbi{a}_\mbiu{E}(i,j) \psi^{l}(i+1,j) + \label{eq:FVM02}\\  
&   & \mbi{a}_\mbiu{W}(i,j) \psi^{l}(i-1,j) + \nonumber\\
&   & \mbi{s}_\mbiu{P}(i,j) \nonumber
\end{eqnarray}
}}

\noindent 
where $l$ represent an iteration. This iteration is incremented when all the lines have been swept. The next algorithm
describes the steps required to solve the system (\ref{eq:FVM02}).

\begin{algorithm}[h!]
\caption{Line $y$ sweep 2D}\label{alg:sweepY2D}
{\footnotesize{
\begin{algorithmic}[1]
  \STATE INPUT: max\_iter, tol, ei, bi, $\psi^0$
  \STATE OUTPUT: Solution $\psi$
  \STATE $l = 0$
  \STATE $a \leftarrow \mbi{a}_\mbiu{P}(i,\_)$; 
  $b \leftarrow \mbi{a}_\mbiu{N}(i,\_)$;
  $c \leftarrow \mbi{a}_\mbiu{S}(i,\_)$; 
  \WHILE{ residual $>$ tol \& $l$ $<$ max\_iter}
  \STATE $d \leftarrow \mbi{a}_\mbiu{E}\psi^{l}(i+1,\_) + \mbi{a}_\mbiu{W}\psi^{l}(i-1,\_) + \mbi{s}_\mbiu{P}(i,\_)$;
  \FOR{$i$ = bi to ei}
  \STATE $\psi^{l+1}(i,\_) \leftarrow $ Thomas($a,b,c,d$)
  \ENDFOR
  \STATE $\psi^{l}(i,\_) \leftarrow \psi^{l+1}(i,\_)$
  \STATE Update residual
  \STATE $l \leftarrow l + 1$
  \ENDWHILE
\end{algorithmic}
}}
\end{algorithm}

This is similar to a block Jacobi iteration, where we take advantage of the tridiagonal
shape of the blocks of (\ref{eq:laplace}) where we apply Thomas algorithm. Note that
the \textbf{for} loop, line 7 in algorithm \ref{alg:sweepY2D}, can be done in parallel
directly.

\section{Numerical Schemes}

 

\end{document}